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The expansions of the spatial distribution of X-rays in suitable systems of polynomials 
appear to converge rapidly. Calculations are performed for various source geometries in an 
infinite homogeneous medium. The surprisingly good convergence of this method of calcu- 
lation is discussed with a view to its possible use in a variety of problems. 



I. Introduction 

Theoretical work at the National Bureau of 
Standards on the propagation of X-rays through 
large masses of homogeneous materials has so far 
been directed primarily to the study of extremely 
deep penetrations [1, 2, 3]. 2 In the complementary 
case of shallow penetration it was thought that a 
straightforward calculation of the effect of successive 
Compton scattering processes would prove con- 
venient. The intention in this paper is to present 
a different method of calculation that can be carried 
out conveniently for moderate penetrations (up to 
10 to 15 mean free paths of the primary radiation) 
and that relies on appropriate application of con- 
ventional physico-mathematical techniques. 

The propagation of X-rays through an infinite 
medium is governed by the transport equation 

u-gradN(r,u,\) = -n(\)N+ f dX'&(\',\) 



i 



du' — 5(1- 

4ir ^7T 



-M-w'-X+XOiV(r,ii',\')+source 

where 

N(r,u,\) is the density of photons of wavelength 
X (in Compton units), at the point r traveling in the 
direction u per unit solid angle; ^i(X) is the total 
narrow beam absorption coefficient; &(X',X) is the 
Klein-Nishina differential coefficient for Compton 
scattering with a wavelength change from X' to X; 
and 8(x) is Dirac's delta-function. 

It is widely known from the study of transport 
equations of type (1) that successive moments of 
the space distribution N are related by a recurrence 
relationship. The zero-th moment, that is, the 
simple space integral of N, has already been calcu- 
lated numerically in a number of sample cases 
[4, 5]. A few more moments can be calculated 
numerically without excessive labor by manual 
operation. Automatic computers might prove use- 
ful for the same purpose.* 



i Work supported by the Office of Naval Research, Mathematics Branch. 
2 Figures in brackets indicate the literature references at the end of this paper. 
* An extensive program of applications is now in progress under the sponsor- 
ship of the Atomic Energy Commission, New York Operations Office. 



It is also widely recognized that a knowledge of 
the moments of a distribution function enables one 
in principle to reconstruct the function itself. How- 
ever, little effort seems to have been applied along 
this line in the study of X-ray or neutron diffusion 
problems, presumably because the method was not 
expected to prove convenient. Actually, the mo- 
ment method can be used very conveniently to 
obtain rapidly convergent expansions of the distri- 
bution function in suitable polynomial systems. 

II. General Discussion of Polynomial 
Expansions 

r 

Consider, for simplicity, a distribution function 
/ of a single space variable x, which ranges from 
— oo to -f oo ; / is assumed to vanish sufficiently 
rapidly at both ends of this range. To obtain a 
polynomial expansion of/, a suitable weight func- 
tion w(x) is first chosen that vanishes at ±oo. 
One then chooses a succession of polynomials 
p , pi(x), . . . p n (x) ... , which has an adjoint 
succession (also polynomials) p+, pt(x) . . ., etc., 
such that 

w(x)p+, (x)p n (x)dx = 5 nn ,. (2) 

J — CO 

One can then write 

j(x)=w(x)2 n a n p n (x)r] (3) 

The coefficients a n are given by the standard formula: 



«„=£ 



pt {x)j{x)dx. 



(4) 



(In connection with each function w(x), one may, 
but need not, use the special set of polynomials p n 
which are identical with their ad joints pi). 

J 'CO 
x m j(x)dx of the 
— CO 

distribution function determines the values of all 
the coefficients a n and thereby the complete ex- 
pansion (3) with the preselected weight function and 
polynomials. If only the first r moments are known, 
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the first r+ 1 coefficients may be determined exact- 
ly, if the (r+ l)-th adjoint polynomial is of degree /. 

[f the function f(x) is the solution of a transport 
equation, it may not be necessary to calculate the 
moments first and thence the coefficients a n . There 
exist recurrence relationships between the successive 
coefficients a n , which are analogous to the relation- 
ships between the successive moments. One may 
then calculate the a n 's one after the other directly. 

The significance of the expansion (3) of a distribu- 
tion function into a series of polynomials is illus- 
trated by the following argument: The first term of 
the series contains a polynomial of grade zero, that 
is, the constant a p . Therefore the weight function 
w{x) itself represents the initial approximation to the 
function /(V) to within the factor a p . If w(x) were 
chosen, just equal to the unknown function f(x) , the 
expansion would reduce to a single term. In gen- 
eral, one can expect the expansion to converge rapidly 
if he chooses a weight function that constitutes a 
sufficiently good approximation to /(#). 

For example, one may think of embodying into 
w(x) advance information, that may be available on 
the behavior of f(x) for very large and very small 
values of x. It is a matter of convenience in any 
particular application how complicated one 4 should 
make w(x) to obtain a rapidly convergent expansion. 
In the practical application to the X-ray problem, 
the choice of the simple exponential function for 
w(x) yields convergence that seems quite adequate 
for a surprisingly large range of x. 

The successive terms of the expansion may be 
described as corrective terms that arc added to the 
initial approximation w(x) «/(#) in order that 
successive moments of the approximation function 

w(x)[a p +a 1 p 1 (x) + a 2 p 2 (x) + ], 



coincide with successive moments of the actual func- 
tion f(x). For example, suppose that an initial ap- 
proximation w{x)~j(x) has been chosen so that 

J 00 /^oo 

w(x)dx= f(x)dx (i. e., a p =l). Then, it 

may turn out that the actual value of the first mo- 
ment I xj{x) dx is larger than the "estimated" 






value xw{x)dx. To correct for this discrepancy 

J -co 

the initial approximation should be modified by the 
addition of a corrective term w(x)a,iPi(x) which in- 
creases w(x) on the side of positive x and reduces it 
on the negative side. Actual calculation of high 
order moments should usually be required only if 
one has to fit the "tails" of the distribution function 
f(x) . It should not introduce substantial corrections 
to the main body of the distribution. 

III. Illustrative Examples 

As a demonstration, we have examined the La- 
guerre polynomial expansions of two otherwise 
known distribution functions, namely exp(— x) 
J (2i-y/x) and exp(— x)iFi(— 3/2, 1,—x), for x>0. 



(Laguerre polynomials are the self-adjoint polyno- 
mials with the weight factor exp (—x) over the range 
to oo). These distributions represent the solutions 
of schematized X-ray or neutron problems [1, 6]; the 
former pertains to the constant mean free path case 
(/x(X) = const.), the second to the case in which 
u(\)=a+b\. 

The values obtained from successive polynomial 
approximations of the "build-up factors" Jo(2i-yJx) 
and [F x {— 3/2, l,—x) are compared with the exact 
values in tables 1 and 2. The fit is surprisingly 
good. 3 The addition of one more polynomial to the 
expansion appears to extend the range of satisfactory 
fit by at least an interval 5x^3 into the tail of the 
distribution function in the direction of increasing x. 
The tables include the results of expansions in 
Laguerre polynomials, as well as expansions in the 
related set of U n polynomials that are introduced 
in the next section. 



Table 1. 



Representation of Jo(2i-y/x) by various finite sums 
of L n (x) and U n (x) polynomials 



X 


Number of Laguerre polynomials 


Exact 


Number of U n 
polynomials 


3 


4 


5 


6 


7 


4 


5 



1 
2 
4 
6 

8 
10 
12 
14 
16 

18 
20 


1.359 
2.039 
4.077 
12. 23 

25.82 

44.85 
69. 32 
99.22 


0. 90(12 
2.341 
4.228 
11.18 
25.37 

50.67 

90. 01 
146. 8 
227.1 
332.0 

465. 7 
632.2 


L.020 

2.270 
4.265 
11.29 

24.81 

49.57 
91.25 

157.8 

261.6 

411.8 

622.8 
910.3 


0. 9968 
2.281 

264.~8~" 
427.9 

669.6 
1018 


L.OOl 

674.3 
1038 


l 

2.280 

4.251 

11.30 

24.88 

49.21 
90. 51 

157.7 

264.0 

427.6 

673.5 
1037 


0. 4797 
2. 545 
4.530 
10.82 
24.47 

50.60 
94. 34 

160. 8 

255. 2 

382.5 

547.9 


0. 9961 
2. 280 
4.289 
11.28 
24.72 

49.27 
91. 47 

159.8 

264. 4 

417.7 

633.7 
928.3 



Table 2. Representation of \F X (—8/2, 1, —a:) by a sum of 
four terms of the Laguerre or U n expansion 



X 


Laguerre 


U n 


Exact 





0. 9953 


0. 9798 


1 


1 


2.681 


2. 688 


2.678 


2 


4. 677 


4.688 


4. 680 


4 


9. 509 


9. 497 


9.512 


6 


15.31 


15.21 


15.29 


8 


21.90 


21.92 


21.89 


10 


29.13 


29. 30 


29.22 


12 


:;.; si 


37. 23 


37.21 


14 


44.72 


45.65 


45.80 


16 


52.73 


54.42 


54.99 


18 


60.66 


63.41 


64. 71 


20 


68.32 


72.49 


74.93 



To illustrate further the workings of these expan- 
sions, one may consider the integral representation of 
the distribution functions as "inverse Laplace 
transforms". In the case of the Bessel function, we 
have 



3 H.Hurwitz, Jr. has informed us that a similar surprisingly accurate fit was 
obtained in neutron calculations following a method of Placzek quoted in 17], p. 
211 and following. 
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±dt 



/-x i rt +i " 

e-*J,{2i^x)=e~ x ^-. 



i /i-ev * 
t 



=e 



ml "''*'•' 



(5) 



where the L n 's are the normalized Laguerre poly- 
nomials. Thus, the expansion in Laguerre poly- 
nomials is obtained by expanding the Laplace 
transform of the distribution function into powers of 
(l—t)/t. This expansion converges rapidly near 
t=l; that is, in the same general region where the 
integrand exhibits a "saddle point" for x^l. Each 
term of the expansion has a pole at the point t=0, 
where the Laplace transform is singular. An expan- 
sion including about four terms affords a good 
approximation of the Laplace transform from t= °° 
(saddle point for x=0) to perhaps t=% (saddle point 
for x=9). 

IV. Calculation of Polynomial Expansions 

The procedure for obtaining a convenient poly- 
nomial expansion of the distribution of X-ray 
photons varies somewhat depending upon the 
geometry of the source. Therefore, a number of 
different cases will be treated separately. 

1. Plane Monodirectional Source: High-Energy 
Photons Only 

The simplest case, perhaps, is that of a parallel 
beam of X-rays of infinite extension that enters a 
material perpendicularly through the plane x— 0, 
The presence of the entrance surface does not cause 
any difficulty in this case because no high energy 
secondary X-rays will be present at the boundary. 

In this case the photon density N (r, u, X) depends 
only on the components x and u x of the vectors r 
and u. The dependence of N upon the direction 
cosine u x may be expressed by means of an expansion 
in Legendre polynomials 

#(*,*», X)S=.^j^ tf,(s,X)P,(tO. (6) 

l 47T 



The incident beam is represented in eq 1 by the 
source term 



Source=g(s) d (1 Ux) /(X), 



(7) 



where /(X) is the incident spectrum. 

Multiplication of eq 1 by P t (u x ), followed by in- 
tegration oyer u yields: 



f X rfX^(X , ,X)P,(l-XfX0^(x,X / ) + 5(x)/(X) (8) 

The deep penetration of X-rays is controlled by 
the absorption of their most penetrating com- 
ponents. [1]. Therefore, one may describe the de- 
pendence of the photon density on the depth of 
penetration x by means of a sum of polynomials in x 
with the weight factor exp (—ax). An obvious 
choice for the parameter a in the weight factor is the 
value of /x(X) for the most penetrating component of 
the radiation, but a less obvious choice of a might 
well prove better. 

As we mentioned earlier, the ortho-normal set of 
polynomials corresponding to the weight factor exp 
(—ax) is the set of Laguerre polynomials given by 
L n (ax)=2 v (n\/(n— v) \v\ 2 ) (—ax)\ We accordingly 
represent the function N t (x, X) in the form 



Ni(z,\y- 



^Zn(X) = 



■r 



>2 n a ln (\)L n (ax), (9) 



L n (ax)Ni(x, \)adx. (9a) 



The equation that determines a Zn (X) is found by 
multiplying eq 8 by aL n (ax) and integrating over x 
from to oo . The integration of the left side of (8) 
is carried out by parts. The values of iV w and Ni-i 
at the limits and <» vanish as shown by the follow- 
ing arguments : 

1. No radiation penetrates to x= oo ; 

2. The density of incoming photons at the entrance 
surface should be taken equal to zero if the source 
term 8 (x) is included in the range of integration; 

3. There are no photons going back out of the 
entrance surface, as the treatment is limited to high 
photon energies. 

The following identity is used in carrying out the 
remaining integration, which involves the derivative 
of the Laguerre polynomial 

ay n '=o 

The integrals of all terms of (8) are thus reduced to 
the form (9a) . The result is j 



^TTT S [(J>+l)a>i+i,n' +la>i-i.n>] = — M(X)a z „(X) + 

r x dx^(x , ,x)P z (i-x+xo^n(x / )+«/(x). (10) 

These equations coincide, for n=0, with those 
whose numerical solution has been discussed in a 
previous paper [5]. The equations for r^>0 can be 
solved in succession by the same method, because 
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Figure 1. Chain linkages of the Legendre-Laguerre coefficients. 

the sum on the left side contains only terms with 

It is often of interest to know the density of pho- 
tons at a point irrespective of their direction. This 
density is represented by the Legendre coefficient 
N (x,\). Calculation of the coefficients a 0n (X) of 
the Laguerre expansion of A 7 requires, according to 
(10), the calculation of a set of coefficients a in with 
Z>0. Figure 1 shows diagrammatically how, for 
instance, the solution of the equation which deter- 
mines a 12 (\) requires a prior knowledge of the func- 
tions a i, a Q0 , a 2 i, and a 20 . In turn, the determina- 
tion of a i and a 2 i requires a prior knowledge of a !0 
and a 30 . If one calculates the moments 



'""Jo 



x n Ni(x,\)dx, 



instead of calculating directly the Laguerre coeffi- 
cients a m, one finds the simpler scheme of interde- 
pendence illustrated by figure 2. Here, solution of 
the equation for b in requires prior knowledge of 
bi-i, n-\ and b i+h n _ij thus, to calculate a 00 , a i, #02, 
and d 3 one must also calculate a !0 , a 2 o, a so> a n, &21, 
and a V2 ) that is, one must solve numerically a total 
of 10 equations This does not involve an excessive 
amount of labor. 

It may be interesting to relate the representation 
of the photon density by means of Laguerre poly- 
nomials to its representation as an inverse Laplace 
transform: 

The moments of the distribution N t in depth cor- 
respond to the coefficients of the Taylor expansion 
of faip, X) in powers of p [7]. The expansion of N t 
in Laguerre polynomials is obtained by expanding 
<t>i into a series of terms of type a n p n l{a—p) nJrl . 



Tl 2 




Ficure 2. Chain linkages of the spatial moments of the Legendre 
coefficients. 
Shaded and unshaded sets of boxes are not interlinked. 



Each of these terms has a pole at p=a. Notice 
that <t>i itself has a singularity at p=Hs, where /x s is 
the smallest value of /*(\) [2]. Thus, the expansion 
of the Laplace transform in powers of p/(a—p) may 
be expected to converge more rapidly than the ex- 
pansion in powers of p if a is given a value near 
that of /x. s . 

2. Plane Isotropic Source 

This case, in which radiation is emitted in all 
directions from all points of a plane, is familiar in the 
study of neutron penetration and diffusion [6, 7]. 

The distribution from a plane source, isotropic as 
well as monodirectional, depends only on the coordi- 
nate x of each point r and on the component u x of the 
direction of flow u. The isotropic source term to be 
entered in eq 1 takes the form: 



Source = 



47T 



(11) 



Correspondingly, one finds, instead of eq 8 

+ rd\'k(\',X)P l (l-y + X)N l (x,\') + 5(x)f(\)5 l0 (12) 
Jo 

If the source is surrounded by an infinitely ex- 
tended medium, the radiation is distributed sym- 
metrically with respect to the source plane x=0 

N(x,u x ,\)=N(-x,—u x ,\); N,(x,\)=(—l)'Ni(-x,\). 

(13) 
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In this case the moments of the space distribution 
involve an integration over the entire range 
— oo <^x<^ oo . Owing to the symmetry conditions 
(13), if Z is even, all the odd moments of Ni are zero 
and if / is odd, all the even moments of N t are zero. 

In the preceding case the nth Laguerre polynomial 
added to the expansion reflected, in essence, informa- 
tion gathered from the calculation of the n-ih moment 
of the distribution function. In the present case 
every other moment vanishes; that is, does not 
contribute to our knowledge of the distribution 
function; therefore, the method applied in the 
preceding case requires some modification. The n-th 
polynominal should now reflect information provided 
by the value of the %n-th or 8n-{-l-th moment. An 
appropriate set of polynomials is required for this 
purpose. 

The weight function should be selected in such a 
way that the distribution function has the correct 
symmetry. To maintain an analogy with the case of 
a monodirectional source, one may choose the weight 
function exp (— a\x\) for the expansion of the 
even-/ Legendre coefficients Ni (x, X). Correspond- 
ingly, for odd values of I the weight function 
xexp(— a\x\) is perhaps the simplest choice. Wewrite 



Ni(x 9 \) = e- a ^^a ln (\)U n (a\x\) 9 I even 



N l (x,\) = axe- a WJla ln {\)V n (a\x\), I odd 

n = 



r (i4) 



where U n and V„ are polynomials of degree n still 
to be determined. 

As indicated in section IV, 1, the first n terms of 
each expansion should embody all the information 
that is provided by the values of the first n usable 
moments of N t (x,\). In other words, the higher 
terms of the expansion (14) should not contribute to 
the lower usable moments. In the case of even I, 
this implies that the component exp (— a\ x |) U n (a\x \) 
of the function Ni (x, X) has vanishing low-order 
even moments 

" x 2n 'e-<* lxl U n (a\x\)dx = for n'<n. (15) 

Similarly: 



/: 



x 2n ' +1 axe~ ah[ V n (a\x\)dx = for n' <n. 



(15a) 



These conditions determine the sets U n and V n to 
within a normalization constant. We find the fol- 
lowing expressions 



Un{y) 2"n\ |_c>3/ + J UMn-m iV 



(16) 



A number of these polynomials are given in table 3 . 4 
Table 3. The first four U n and V„ polynomials 



U = l 


v| 


u i=\d-y) 


Vi-\@-y) 


^2=| (3— 52/-{-2/ 2 ) 


^2=^ (15— 92/ + 2/ 2 ) 


u *=k (M-^y+^-y 3 ) 


F 3=3^ (t05-872/-hl82/ 2 -2/ 3 ) 



To determine coefficients of the expansion (14) 
one may utilize the sets of adjoint polynomials UZ> 
and y+/ , which fulfill the orthogonality condition (2) 
together with U n and V n , respectively. These 
polynomials are 



nW hi (2,-1)! V v ) y J 



(17) 



We multiply eq 12 by aUt (a x), if I is even, or by 
aV+ (ax), if I is odd, and integrate over x from 
— oo to oo . The integration proceeds as the integra- 
tion of (8) . The left side integrated by parts yields 
one term that vanishes and an integral that contains 
the derivative dU+/dx (I even) or dV+/dx (I odd). 
These derivatives fulfill the identities 



dy V "- lW 

dVi(y) 



dy =S^. 



(18) 



Thereby, the integration of eq 12 reduces to inte- 
grals of the type 



am(X)= \ Ui{ccx)Ni{x,\)dx 

= §C)w J "" Wi 

a,„(X)= f" V:(ax)N,(x,\)dx 



I even 



-^Y'H-i y-ir 1 h en 
~h\ v )j2^i)i" 2 '- i{x) ' 



I odd 



(19) 



4 It might be worth mentioning in passing that U n and V n are solutions of 
the following 3d order differential equation, with the parameter m equal to zero 
and one, respectively 

d*W d*W dW 

^^ + (2ra+1 ~ 3l/) l^ +2(?/ - 2w - 1) ^- 27lTy=0 

Higher associated polynomials exist having weight functions and adjoint poly- 
nomials analogous to those of U n and V n . 



450 




Figure 3. — Interlinkages for isotropic sources. 
Shaded elements only are involved in the calculation. 

where bi tn f is again the 7i/-th spatial moment of 

Ni(x,\). ' The result is 



f X da; , fc(X',X)Pz(l--X+X0az w (X0 + a/(X)6 i o, Z 
Tw-rT S [(J+l)a,+i, n' + ^i_i t «'] = — M(X)a^(X) + 

^6-j- 1 n' = 

f dX'ifX',X)P,(l--X+X')ai ll (X'), 'odd 

(20) 



even 



These equations can be solved chainwise, much like 
the eq 10. Notice that a tn = for n<(l/2 — l), 
since the source term is confined to the equation 
with Z=0. The equations must be solved in the 
following order: a m a 10 , a 2ll a 3h • • •, a i,On, a 2i • • •, 
(Lo2j «i2. • • •» &03 • • •- Here again, as in the case 
of the moaodirectional source, the calculation of four 
terms of the expansion of iVo (x, X) requires the 
solution of 10 equations. 

The fact that a Zw =0 forn< (1/2 — 1) relates to the 
known fact that the moment 6,^= I x n Ni(x, X)dx 

J— CO 

vanishes when n<Cl [7]. (Only the moments corre- 
sponding to the shaded boxes in fig. 3 are not zero. 
The value of each moment bi n depends on the val- 
ues of the moments b t 



+ l,n-l> 



and b? 



as shown 



in fig. 2. The calculation must begin with & o, 
whose equation is inhomogeneous, that is from the 
upper left corner of the scheme.) Accordingly, the 
Laplace transform of Ni(x,\), namely, <l>i(p,\), has 
a Taylor expansion that begins with the Z-th term. 
The polynomial expansions (14) of Ni are obtained 



as inverse Laplace transforms if <f> is represented as 
a series of terms p 2n /(a 2 —p 2 ) n + l , for even I, or 
p 2n + l l(ol 2 -p 2 ) n + l for odd I 

3. Plane Monodirectional Source No Energy 
Limitation 

The method developed for the treatment of X-rays 
from a plane isotropic source can be applied to the 
study of radiation from other plane sources. To 
this end, one must generalize the preceding treat- 
ment somewhat. 

It will be remembered that because the source is 
symmetric and because of the interlocking of mo- 
ments illustrated in figure 2, the only useful moments 
are those which correspond to the black squares of 
the checkerboard pattern of figure 2. The half- 
checkerboard pattern of useful moments (fig. 3} 
results from the isotropy of the source (eq 11). In 
the general symmetric case, source terms, will be 
present in all equations, not just those alone in 
which /=0; and the full checkerboard scheme of 
moments (fig. 2) will be used. 

On the other hand, if the source were antisym- 
metric in x, the moments corresponding to the 
shaded squares in figure 2 would all be zero. The 
useful moments would be interrelated among them- 
selves just as before and would correspond to the 
white squares of figure 2. 

One can take advantage of this separate inter- 
locking of different sets of moments by regarding a 
plane monodirectional source as the superposition 
of a symmetrical and an antisymmetrical source. 
The symmetrical source component injects equal 
photon densities in opposite directions perpendicular 
to the plane x = on either side. The antisym- 
metrical source injects a photon density equal to 
that of the symmetrical source of the positive side 
of x=0 and a "negative" photon density of equal 
size on the negative side. 

The photon distribution N {a) (x,u x ,\) due to the 
antisymmetric source has symmetry properties 
opposite to (14), namely, 

N^(x,u x ,\)=-N ia) (-x,-u x ,\); 

N l ia) tx,\) = (-l) l+l N l (a) {-x,\). (21) 

The polynomial expansions of the Legendre coeffi- 
cients iV7 o) (x, X) involve the V n polynomials for even 
I and the U n polynomials for odd I, that is, they 
conform to (14) with opposite even-odd condition. 
The total photon density from a plane mono- 
directional source is thus given by the sum of two 
terms 

N(;x,u x ,\)=N^(;x,u x ,\)+N^(x,u x ,\). (22) 

N (s) is determined by a set of equations identical with 
(20), except that the source term does not include 
the factor 8 i0 . N {a) is determined by a similar 
system with even and odd I interchanged. 

The total density N(x, u x , X) all but cancels on 
the negative side of the source plane, inasmuch as 
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N (a) is negative there. The residual photon density 
for z<0 represents radiation that has been back- 
scattered behind the source. It should consist of 
comparatively low energy photons. 

4. Point Isotropic Source 

The photon density generated by a point isotropic 
source can be derived from the solution corresponding 
to a plane isotropic source (case (b)) [7]. This 
derivation is especially simple if one calculates only 
the total density of photons at each point irrespective 
of their direction. For this reason we shall outline 
here only sketchily an independent calculation for 
the case of a point source. 

Owing to the spherical symmetry of the problem, 
the distribution function N(r, u, X) depends only on 
the magnitude r of r and on the component u r =u-r/r 
of u, if the origin of the coordinates is laid at the 
source. The source term in eq 1 takes the form 
8 (r)/(X)/4 7rr 2 and the transport eq 1 becomes 



u • gradN(r, u T ,\)=u r 



Z>N . (1- 



dr 



~du T 



i 



■ix(\)N+ rd\'k(\',\) 



du' ^8(1-11- u'-\ + y) 
iv{r,u r ,\)-t 4 ^ 2 (23) 

The density N can be expanded in Legendre poly- 
nomials of u r , as in eq 6. The transport equation 
takes then a form corresponding to (8), whose deri- 
vation makes use of the identity 

a-u*)dp l /du r =i(i+i)(p l - l -p l+l )/(2i+i). 

=-»(\)N l (r > \)+rd\'k(X,*.)P l (l--k+\')N l (r,\') 

■ S(r)/(\)fto 
"•" 4xr 2 (24) 

In this case we do not derive directly the coeffi- 
cients of an expansion of Ni(r, A) in polynomials in r, 
but we calculate the moments 



I 



r n Ni(r,\)4:wr 2 dr=bi n (\). 



To this end one multiplies (24) by r n and integrates 
over all the space. The result is 



2Z+1 



[(l+m-n)b l+l , n - 1 (\)-l(l+n+l)b l . ltn ^(\)] = 



-m(X)MX) + 
+/(X)5jo<5 w o 



j: 



d\'fc(X',X)P,(l-X + X')&in(A') 



(25) 



This set of equations can be solved chainwise in 
exactly the same way as in the plane isotropic case. 
Coefficients of the U n and V n polynomials can be 
obtained by adding and subtracting moments. The 
coefficients in the expansion of the total density of 
photons in U n polynomials are related in the follow- 
ing way with the corresponding coefficients in the 
plane isotropic source case 

(tto»)poliit = (2n-{- l)(& rc)piane ~~ ^ n ( a 0, n - l)plane- (2 6) 

The coefficients a Xn for the two sources in the case 
1^0 are related in a more complicated way. 



5. Point Monodirectional Source 

This kind of source may properly be regarded as 
the elementary source, since any source can be 
represented as an aggregate of point monodirec- 
tional sources. The case of a point monodirectional 
source is more complicated than the cases considered 
above, because it involves explicitly the diffusion of 
radiation sidewise from the initial "line of fire." 

Cylindrical coordinates r=(z, p, </>), with the z 
axis lying along the line of fire are appropriate to the 
geometry of this problem. The direction vector u 
can be described by polar coordinates (0, $+<£), so 
that the zero azimuth $—0 refers to a plane through 
the z axis and the point r under consideration. 
Owing to the cylindrical symmetry of the source the 
photon density depends explicitly on the variables 
in the following way 



N(r,u,\)=N(z, P ,e,*,\). 



(27) 



To calculate the expression u- grad N explicitly, 
one must consider that the gradient is taken for 
constant u, which means 8$=— 8$. Accordingly 

ugrad N=cos 0^- + sin cos <S>-s— 

OZ Op 



sin 6 sin $ dN 
P d$ 



(28) 



In this problem the distribution function should be 
expanded into Laplace spherical harmonics F* m (0, <£) 
and moments should be taken with respect to both 
space coordinates z and p. Therefore we introduce 
the set of dependent variables 

N?£(\) = fV+Wp [ " z n dz f'sinfldfl 
Jo J-» Jo 



f 



d*Y? m (0,*)mz, P ,e,$,\) 



(29) 



A corresponding four-fold integration applied to eq 
1 reduces this equation to the form 
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■-= =\n^(l+l) 2 -m 2 Nj+t£ 

V4(/+l) 2 -l( 

-2±^ V(i-m+l)(Z-m + 2) N^-i 
+ lZ^L V(Z+m+l)(/ + m + 2) tfftf r.Vi 



V4/ 2 -l 



n^/P-m 2 Nllti 



+ 2 ^ ja+rn)Q+m-l)NTLrt-i 



p — m 



V(f-™)(J-m-l)JVriVm+i 



=- M (x)i^r^(x)+J^dx / k(x / l x)P I (i-x+xoJVw(x / ) 

+/(X)5„ 5 P o5mo -\~^r (3 ) 

This equation reduces to the equation for the 
plane monodirectional source, as it should do, when 
p = m = 0. The coefficients JVf t b°(X) can be calcu- 
lated accordingly. Once this is done one can pro- 
ceed to p = l, m=±l, and then on again. 

In the preceding cases the calculation of a three- 
variable function such as N (r, w r , X) was reduced to 
the calculation of a two-dimensional array of co- 
efficients Ni(\). The present problem involves 
five independent variables and reduces to the solu- 
tion of an "array of problems," each of them as 
complicated as any of the previous ones and corres- 
ponding to one pair of values of p and m. The 
genetic connections in the array of problems are 
represented in figure 4. 

Various polynomial representations of the radial 
distribution might suggest themselves in the course 
of any practical application. The use of a gaussian 
weight function should presumably be convenient 
at least in the lower energy range. 

V. Numerical Applications* 

1. Penetration of 10-Mev Gamma Rays in Lead 

The method of section IV, 1 was applied numeri- 
cally to the case of a parallel beam of 10.2-Mev 
photons entering an infinitely thick layer of lead. A 
heavy element was chosen because it is a more diffi- 
cult situation to treat with asymptotic methods, 
owing to the minimum of the absorption coefficient 
near 3 Mev [2]. The same initial energy was 
chosen which is used by Karr and Hurwitz in their 
asymptotic calculation with the "straight ahead" 
approximation [8]. The calculation was extended 
only down to photon energies of about 1 Mev 
(=0.49), partly because of the limitation to the 
validity of the method indicated in section IV, 1. 



* Additional applications were reported by L. V. Spencer and F. A. Stinson 
at the New York Meeting of the American Physical Society, February 3, 1951. 



This calculation was carried out not omy by the 
method of section IV, 1 but also with a correspond- 
ing "straight ahead" approximation in order to esti- 
mate the error involved in neglecting angular 
effects. The straight ahead approximation implies 
disregarding the effects of deflections experienced by 
high-energy photons in the course of Compton scat- 
tering. In this approximation one solves only the 
eq 8 for 1=0, replacing diVi/dz by &N /c)x, so that 
iVo (x, X) is the only distribution function to be 
determined. 

Finally, the straight ahead calculation was car- 
ried out once more, neglecting the sin 2 terms in 
the Compton cross section, that is, omitting the 
last two terms of the expression 



w.«-!(x)"[y 



+ y~ 2(X — X') + (X — X') 2 JMThomsol. 



2 ]/ 



The purpose of this calculation was partly to esti- 
mate the importance of these terms and partly to 
obtain a proper comparison with the work of Karr 
and Hurwitz, who used this approximation in mak- 
ing their calculation. 

Four terms of the polynomial representation (9) 
of the function N (x, X) were worked out with a 
set equal to fj, my the minimum value of /*(X). For 
comparison, the coefficients were regrouped to cor- 
respond to a polynomial expansion with a set equal 
to ju(X=0.05) and also to a polynomial expansion 
with the weight factor x exp (— /j, m x), the last being 
perhaps the most realistic weight factor in view of 
what is known about the initial and asymptotic 
behavior of the build-up factor. The results using 
these different weight factors were all in satisfactory 
agreement. 

Figure 5 shows the differential spectrum given by 
calculation according to section IV, 1. The energy 
density is plotted as a function of the energy for 
various penetrations. Notice the accumulation of 
energy in the ranges 3 Mev, corresponding to the 
minimum absorption coefficient. 



P 3 




Figure 4. Interlink ages among the radial moments and the 
Fourier coefficients of the azimuthal distribution for a point 
monodirectional source. 

Only shaded elements are involved in the calculation. 



453 



1.000 — 




10 E 



Figure 5. Differential spectra I(E, /j. m x) of the X-ray intensity 
at various distances from a plane monodirectional 10.2-Mev 
source in Pb. 

The position of the source spectrum is indicated by the solid vertical line, its 
intensity by the area of the dotted rectangle. The scale of ordinates is normalized 
to unit strength of the source atz=0; at greater depths it discounts an exponential 
decay corresponding to the absorption coefficient n (3.2 Me v)=/i W = 0.469 cm- 1 of 
the most penetrating component. The dotted line for x=0 departs from the base 
line because only a finite number of terms in the expansion were calculated. 
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Figure 7. Differential spectra I{E, n^r) of the X.-ray intensity 
at various distances r from a Co 60 point isotropic source in 
water. 

The positions of the two lines of the source spectrum are indicated by the solid 
vertical lines, their intensities by the area of the dotted rectangles. The scales of 
ordinates are normalized to unit intensity of the 1.33-Mev component of the 
source; /i =0.0612 cm -1 is the narrow-beam absorption coefficient of this com- 
ponent. 



Figure 6. Differential intensity l(E,fi m x) of the 1.04-Mev 
component of the spectrum of figure 5 as a function of the 
distance from the source, calculated by different methods: 
(a) as in figure 5 ("exact"), (b) disregarding the deflection of 
scattered photons ("straight ahead"), (c) using the extreme 
relativistic form of the Klein-Nishina cross section (also 
"straight ahead"). 

The point at the end of the top curve is the result of the Karr-Hurwitz 
asymptotic calculation. 
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Figure 8. Semilogarithmic plot of the spectra of figure 7, 
which illustrates the trend toward equilibrium. 
The higher energy portions of the spectra have been omitted for clarity. 

Figure 6 shows the angular calculation compared 
with the two straight ahead calculations. The energy 
density of the component A = 0.49 is plotted against 
depth of penetration. Karr and Hurwitz carried 
their asympotitic calculations down to fx m x=lO. It 
is possible to compare their results with ours at that 
penetration. The circle gives their value and is to 
be compared with our upper curve. The two values 
are much closer together than was expected from the 
accuracy of the calculations. Notice that the angular 
calculation differs by a factor of almost 3 from the 
Karr-Hurwitz calculation at 10 mean free paths of 
the most penetrating component. 

2. Co 60 Isotropic Point Source in Water 

This problem has been investigated experimentally 
by White [9]. We have made numerical calcula- 
tions by using the method outlined in section IV, 4. 
The coefficients of the first 4 U n polynomials in 
the expansion of N Q (r,\) have been evaluated. 
(Below 0.150 Mev, the coefficients of U 2 and U 3 were 
estimated from the trend of the lower coefficients, 
since approximate equilibrium was attained.) 

Some results are illustrated in figures 7, 8, and 9. 
Figure 7 shows the differential energy spectrum at 




5 10 15 

DISTANCE FROM SOURCE , u Q r 

Figuee 9. Comparison of the calculated curve for the build-up 
factor with the experimental results of [9]. 

various depths of penetration. The same results are 
plotted scmilogaritlimically in figure 8 to illustrate 
the approach to an equilibrium spectrum at the 
lower energies. At equilibrium the curves become 
parallel, indicating that the spectrum no longer 
changes with increasing penetration. 

Figure 9 compares White's experimental build- 
up factors, as measured by an ionization cham- 
ber (and a counter), with a theoretical build-up 
factor. The latter was evaluated by weighting the 
calculated spectra with the theoretical response in 
roentgen of a standard free-air ionization chamber 
and then again with the specific correction factor of 
the Keleket III pocket dosimeter [10]. The integral 
response over all spectral components was divided 
by the calculated response to the primary radiation 
to give the theoretical build-up factor. The slight 
downward curvature at the end of the theoretical 
curve is an artifact that indicates that the 4-term 
approximation is beginning to fail at this distance. 

VI. Discussion 

The analytical developments and the numerical 
applications presented in this paper seem to give 
some confidence in the following conclusions: 

a. A method is available for calculating the dis- 
tribution of X-rays in uniform media with a moderate 
amount of labor up to fairly large depths of penetra- 
tion. 

b. This technique can presumably be applied with 
success to the broad class of transport phenomena 
governed by Boltzmann-type equations. It can 
take advantage even of rather crude theoretical 
predictions on the behavior of a distribution function 
to formulate an initial approximation that can 
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thereafter be rapidly improved by straightforward 
numerical work. 

Nevertheless, much further practical experience 
may be required before one can assess with some 
confidence the criteria for the most economical 
application of this technique. For example, it is not 
immediately apparent how far one can conveniently 
push the search for an improved weight function. 

The present method should probably prove con- 
venient for mapping out the distribution of various 
radiations within thick layers of scattering and 
absorbing materials. The necessary formulae for 
calculating the spacial moments of electron distribu- 
tions have recently been developed by Lewis [11]. 

One can conceive of making still broader use of the 
methods considered in section II. There is some 
question as to whether the fitting of distribution 
functions by rapidly converging approximations 
should remain a problem by problem affair or 
whether some useful guiding principles can be 
formulated. 

The present approach is complementary to the 
earlier papers [2, 3] which dealt primarily with the 
asymptotic behavior of the photon distribution and 
which made extensive use of the method of Laplace 
transforms. The work reported in those papers has 
now been developed into a more comprehensive 



treatment which will be the object of a separate 
report. 

Continuation of this work will aim primarily in two 
directions : 

a. To evaluate the effect of boundaries between 
two media and, in general, to consider the propaga- 
tion through inhomogeneous barriers. 

b. To include the effect of X-ray regeneration by 
the cascade shower mechanism, which becomes 
important at high energies. 
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